#ifndef _tppmc_Particle_H_
#define _tppmc_Particle_H_

//
// Particle.h
//
// Created by TB on 11/17/11.
// Copyright 2011 BNL. All rights reserved.
//

#include <cmath>
#include <cstdlib> // For abs(int)
#include <iostream>
#include <vector>

#include <TObject.h>
#include <TVector3.h>
#include <TParticle.h>

namespace tppmc {
   
  class Particle : public TObject {
  public:
    virtual TVector3 PxPyPz() const = 0;
    virtual int Type() const = 0;
    ClassDef(tppmc::Particle, 1)
      };


  /**
     Parton.
     \remark Formerly known as pythpart.
     \todo Add energy field
  */
  class Parton : public Particle {

  public:

    /** Constructor */
    Parton(double momX = 0., double momY = 0., double momZ = 0.,
	   int species = 0, int nucleon = 0, int partonspin = 0,
	   int protonspin = 0);

    /** Destructor */
    virtual ~Parton() { }

    TVector3 PxPyPz() const;

    void SetPxPyPz(const TVector3&);

    /**
       Whether this is the first or second parton
       in the hard interaction.
    */
    int Group() const;
      
    /** The PDG code of the parton */
    int Pdg() const;
      
    /** Particle ID according to LHAPDF convention */
    int Lhapdf() const;
      
    bool IsGluon() const;

    bool IsQuark() const;
      
    /** x component of momentum */
    double Px() const;
    /** y component of momentum */
    double Py() const;
    /** z component of momentum */
    double Pz() const;
    /** Transverse component of momentum */
    double Pt() const;

    /** Longitudinal momentum fraction */
    double X() const;
      
    virtual operator int() const;
      
    /** Parton type */
    int Type() const;
      
    /** Magnitude of 3-momentum */
    double P() const;
      
    /** Parton spin direction */
    int Spin() const;

    /** Parent proton spin direction */
    int ProtonSpin() const;

    /** Sivers weight = sivers / unpolarised PDF */
    double SiversWeight() const;

    /** Sivers spin factor = proton spin * parton px / parton pT */
    double SiversSpinFactor() const;

    void Speak(std::ostream& os) const;
      
    int type;      ///< Parton species
    int group;     ///< from +z nucleon = 1, -z nucleon = 2
    int spin;      ///< spin orientation of parton
    int pspin;     ///< spin of proton
    double px;     ///< x component of momentum
    double py;     ///< y component of momentum
    double pz;     ///< z component of momentum
    double x;      ///< Momentum fraction
    double cos_spin_depol; // cosine of angle between spin axis in parton frame and parton momentum
    double phi_spin; // Azimuthal angle of the spin axis in this parton's frame
    double mF1;
    double mSivers;
    double mH1;
    double mSiversWeight;
    double mSiversSpinFactor;

    ClassDef(tppmc::Parton, 1)
      };
   
  typedef Parton pythpart; // Temporary for backward compatilibility

  inline TVector3 Parton::PxPyPz() const {
    return TVector3(Px(), Py(), Pz());
  }

  inline void Parton::SetPxPyPz(const TVector3& v) {
    px = v.Px();
    py = v.Py();
    pz = v.Pz();
  }

  inline void Parton::Speak(std::ostream& os) const {
    os << "Parton " << type << std::endl;
  }

  inline int Parton::Group() const { return group; }

  inline int Parton::Pdg() const { return type; }

  inline bool Parton::IsGluon() const { return Pdg() == 22; }

  inline bool Parton::IsQuark() const {
    return ::abs(Pdg()) < 7 and Pdg() not_eq 0;
  }

  inline double Parton::Px() const { return px; }

  inline double Parton::Py() const { return py; }

  inline double Parton::Pz() const { return pz; }

  inline double Parton::Pt() const {
    return sqrt(pow(Px(), 2.) + pow(Py(), 2.));
  }

  inline double Parton::X() const { return x; }

  inline Parton::operator int() const { return type; }

  inline int Parton::Type() const { return type; }

  inline double Parton::P() const {
    return sqrt(pow(Px(), 2.) + pow(Py(), 2.) + pow(Pz(), 2.));
  }

  inline int Parton::Spin() const { return spin; }

  inline int Parton::ProtonSpin() const { return pspin; }

  inline double Parton::SiversWeight() const { return mSiversWeight; }

  inline double Parton::SiversSpinFactor() const { return mSiversSpinFactor; }

  /**
     A particle generated from a fragmenting parton.
     \remark Rename class to Hadron
     \todo Write out PDG code of hadron
     \todo add energy field
     \todo fields like jx, jy, jz, group, parton etc are now superfluos with the
     mParton field, so remove them
     \todo remove "accepted_..." fields (PHENIX-specific)
  */
  class Hadron : public Particle {
      
  public:
      
    Hadron();
      
    virtual ~Hadron() { }
      
    void Speak(std::ostream& os) const {
      os << px << " " << py << " " << pz << std::endl;
    }
      
    virtual int Type() const { return type; }
      
    virtual operator int() const { return type; }
      
    virtual double Pt() const { return ::sqrt(::pow(Px(), 2.) +
					      ::pow(Py(), 2.)); }
      
    virtual Hadron& operator+=(const Hadron& other) {
      px += other.px;
      py += other.py;
      pz += other.pz;
      z  += other.z;
         
      return *this;
    }
      
    /**
       Whether the hadron came from the first or second parton
       in the hard interaction.
    */
    int Group() const { return mParton.group; }
      
    /** x component of momentum */
    double Px() const { return px; }
    /** y component of momentum */
    double Py() const { return py; }
      
    /** z component of momentum */
    double Pz() const { return pz; }
      
    /** The momentum three-vector */
    TVector3 PxPyPz() const { return TVector3(Px(), Py(), Pz()); }

    /** Momentum fraction carried by the hadron */
    double Z() const { return z; }
      
    /**
       Generated phi of the hadron i.e. wrt parton spin in the parton frame
    */
    double PhiGen() const { return phigen; }
      
    /**
       Returns true if this was the last, momentum-conserving, hadron in the
       event, i.e. was not generated using fragmentation functions.
    */
    bool IsLast() const { return lasthad not_eq 0; }
      
    /** Returns the Parton from which this Hadron was produced */
    const Parton& GetParton() const { return mParton; }
      
    /**
       Returns the value of the fragmentation function used in generating
       this hadron, including any Collins amplitude.
    */
    double Frag() const { return mFrag; }
      
    /**
       Returns the value of the Collins amplitude for this event.
    */
    double Collins() const { return mCollins; }
      
    /** Returns the Collins angle of this hadron */
    double CollinsAngle() const { return collinsangle; }
      
    int type; ///< Hadron species
    int lasthad; ///< 1 if it was generated to conserve momentum, 0 otherwise
    int favoured; ///< Needed?
    double px; ///< x component of momentum
    double py; ///< y component of momentum
    double pz; ///< z component of momentum
    double z; ///< z of hadron compared to parton
    double collinsangle;
    double phigen;
    double mFrag;
    double mCollins;
    Parton mParton; ///< The parton from which the hadron was produced
      
    static Hadron Test();
      
    ClassDef(tppmc::Hadron, 1)
      };


  /**
   */
  class softpart {
      
  public:
      
    double px; ///< x component of momentum
    double py; ///< y component of momentum
    double pz; ///< z component of momentum
    int type; ///< Parton species
    double e;
  };


  /**
   */
  class jetstruct {
      
  public:
      
    float thrust; 
    double px; ///< x component of momentum
    double py; ///< y component of momentum
    double pz; ///< z component of momentum
    std::vector<Hadron> constituents;
    std::vector<int> hindex;  
  };
   
   
} // namespace tppmc

#endif
